/** 
 * \file gdrotmg.cu
 * \author Kyle E. Niemeyer
 * \date 10/04/2011
 *
 * Based on "drotmg.f" from BLAS.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gdrotmg computes elements of a modified Givens plane rotation matrix.
 * Uses unrolled loops for increments equal to one.
 *
 * \param[in,out] dd1      first diagonal element
 * \param[in,out] dd2      second diagonal element
 * \param[in,out] dx1      x-coordinate
 * \param[in]     dy1      y-coordinate
 * \param[out]    dparam  rotation matrix information, dimension 5
 */
__device__ __inline__ void gdrotmg ( double *dd1, double *dd2, double *dx1, double dy1, double *dparam )
{
  
  const double gam = 4096.0;
  const double gamsq = 16777216.0;
  const double rgamsq = 5.9604645e-8;
  
  double dflag, dh11, dh12, dh21, dh22;
  
  if ( dd1 < 0.0 ) {
    
    dflag = -1.0;
    dh11 = 0.0;
    dh12 = 0.0;
    dh21 = 0.0;
    dh22 = 0.0;
    
    dd1 = 0.0;
    dd2 = 0.0;
    dx1 = 0.0;
    
  } else {
    
    double dp2 = dd2 * dy1;
    
    if ( dp2 == 0.0 ) {
      dflag = -2.0;
      dparam[0] = dflag;
      return;
    }
    
    double dp1 = dd1 * dx1;
    double dq2 = dp2 * dy1;
    double dq1 = dp1 * dx1;
    
    if ( fabs(dq1) > fabs(dq2) ) {
      
      dh21 = -dy1 / dx1;
      dh12 = dp2 / dp1;
      
      double du = 1.0 - dh12 * dh21;
      
      if ( du > 0.0 ) {
        dflag = 0.0;
        dd1 /= du;
        dd2 /= du;
        dx1 *= du;
      }
      
    } else {
      
      if ( dq2 < 0.0 ) {
        dflag = -1.0;
        dh11 = 0.0;
        dh12 = 0.0;
        dh21 = 0.0;
        dh22 = 0.0;
        
        dd1 = 0.0;
        dd2 = 0.0;
        dx1 = 0.0;
      } else {
        dflag = 1.0;
        dh11 = dp1 / dp2;
        dh22 = dx1 / dy1;
        du = 1.0 + dh11 * dh22;
        double dtemp = dd2 / du;
        dd2 = dd1 / du;
        dd1 = dtemp;
        dx1 = dy1 * du;
      }
    }
    
    if ( dd1 != 0.0 ) {
      while ( ( dd1 <= rgamsq ) || ( dd1 >= gamsq ) ) {
        
        if ( dflag == 0.0 ) {
          dh11 = 1.0;
          dh22 = 1.0;
          dflag = -1.0;
        } else {
          dh21 = -1.0;
          dh12 = 1.0;
          dflag = -1.0;
        }
        
        if ( dd1 <= rgamsq ) {
          dd1 *= gam * gam;
          dx1 /= gam;
          dh11 /= gam;
          dh12 /= gam; 
        } else {
          dd1 /= gam * gam;
          dx1 *= gam;
          dh11 *= gam;
          dh12 *= gam;
        }
      }
     
    }
    
    if ( dd2 != 0.0 ) {
      while ( ( fabs(dd2) <= rgamsq ) || ( fabs(dd2) >= gamsq ) ) {
        
        if ( dflag == 0.0 ) {
          dh11 = 1.0;
          dh22 = 1.0;
          dflag = -1.0;
        } else {
          dh21 = -1.0;
          dh12 = 1.0;
          dflag = -1.0;
        }
        
        if ( fabs(dd2) <= rgamsq ) {
          dd2 *= gam * gam;
          dh21 /= gam;
          dh22 /= gam;
        } else {
          dd2 /= gam * gam;
          dh21 *= gam;
          dh22 *= gam;
        }
        
      }
    }
    
  }
  
  if ( dflag < 0.0 ) {
    dparam[1] = dh11;
    dparam[2] = dh21;
    dparam[3] = dh12;
    dparam[4] = dh22;
  } else if ( dflag == 0.0) {
    dparam[2] = dh21;
    dparam[3] = dh12;
  } else {
    dparam[1] = dh11;
    dparam[4] = dh22;
  }
  
  dparam[0] = dflag;
  
}